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Abstract 

In the paper, we propose an algorithm for the efficient parallel implementation of clastoplastic prob- 
lems with hardening based on the so-called TFETI (Total Finite Element Tearing and Interconnecting) 
If) domain decomposition method. We consider an associated elastoplastic model with the von Mises plastic 

criterion and the linear isotropic hardening law. Such a model is discretized by the implicit Euler method 
in time and the consequent one time step elastoplastic problem by the finite element method in space. 
^ The latter results in a system of nonlinear equations with a strongly semismooth and strongly monotone 

y-— operator. The semismooth Newton method is applied to solve this nonlinear system. Corresponding 

linearized problems arising in the Newton iterations are solved in parallel by the above mentioned TFETI 
domain decomposition method. The proposed TFETI based algorithm was implemented in Matlab par- 
allel environment and its performance was illustrated on a 3D elastoplastic benchmark. Numerical results 
for different time discretizations and mesh levels are presented and discussed and a local quadratic con- 
vergence of the semismooth Newton method is observed. 

<N 

. £r 1 Introduction 

Q\ Elastoplastic processes describe behaviour of solid continuum beyond reversible elastic deformations. They 
are typically described by hysteresis models with a time memory [BS961 IKr96| . The rigorous mathematical 
analysis of elastoplastic problems and the numerical methods for their solution started to appear in the late 

^ 70ies and in the early 80ies by the work of C. Johnson [Jo76l lJo78] , H. Matthies |Ma79l IMa79b| , V. Korneev 
and U. Langer [KL84J, J. Necas and I. Hlavacek |NH8L and others. Since then a lot of mathematical 
contributions to computational plasticity have been written. Here we refer at least to some of them: to the 

•'-j monographs by J. Simo and T. Hughes [SH98] and W. Han and B. Reddy |HR99j . to the book of Blaheta 

rjj |B199j . to the habilitation theses by C. Carstensen |C93| and C. Wieners |Wi00] . and to the collection of 
c3 papers by E. Stein et al. [St03j. 

In this paper, we focus on the efficient parallel implementation of elastoplastic problems based on the 
TFETI domain decomposition method. More specifically, we consider an associated elastoplasticity with 
the von Mises plastic criterion and the linear isotropic hardening law (see e.g. [HR99, B199, NPO08]). 
The corresponding elastoplastic constitutive model is discretized by the implicit Euler method in time and 
consequently a nonlinear stress-strain relation is implemented by the return mapping concept (see e.g. 
|NPO081 IACZ 99, B199]). This approach together with the balance equation, the small strain assumption 
and a combination of the Dirichlet and Neumann boundary conditions leads to the solution of a nonlinear 
variational equation with respect to the primal unknown displacement in each time step. Such an equation 
can also be equivalently formulated as a minimization problem with a potential energy functional (see e.g. 
[GV091 [Sy09p . 

* corresponding author, email: jan.valdman@vsb.cz 
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By a finite element space discretization of the one time step problem, we obtain a system of nonlinear 
equations. The corresponding nonlinear operator is nondifferentiable but strongly semismooth. Therefore, 
it is suitable to choose the semismooth Newton method for solving the system since the strong semismooth- 
ness together with other properties ensure local quadratic convergence. Semismooth functions in finite 
dimensional spaces and the semismooth Newton method were introduced in |QS93|. In elastoplasticity, the 
semismoothness was investigated for example in |GV09j ISW11| |Sy09[ |Syl2| . 

In each Newton iteration, it is necessary to solve the respective linearized problem. Different linear 
solvers including those based on multigrid have been successfully tested in |WilO|, IGKLSVllj . Moreover, 
since the linear systems of equations corresponding to the elastic and elastoplastic problems are spectrally 
equivalent |KLV04j . all preconditioners for elastic problems can be applied to elastoplastic ones as well. 

A linear solver considered in this paper is based on a FETI type domain decomposition method enabling 
its efficient parallel implementation. The standard FETI method (FETI-1) was originally introduced by 
Farhat and Roux [FR92]. Using this approach, a body is partitioned into non-overlapping subdomains, an 
elliptic problem with Neumann boundary conditions is defined for each subdomain, and intersubdomain field 
continuity is enforced via Lagrange multipliers. The Lagrange multipliers are evaluated by solving a relatively 
well conditioned dual problem of a small size that may be efficiently solved by a suitable variant of the 
conjugate gradient algorithm. The first practical implementations exploited only the favorable distribution 
of the spectrum of the matrix of the smaller problem |Ro92| . known also as the dual Schur complement 
matrix, but such algorithm was efficient only with a small number of subdomains. Later, Farhat, Mandel, 
and Roux introduced a "natural coarse problem" whose solution was implemented by auxiliary projectors 
so that the resulting algorithm became in a sense optimal I 'M HO 1. IRF98j . Here, we use the Total-FETI 
(TFETI) [DHK06J variant of FETI domain decomposition method, where even the Dirichlet boundary 
conditions are enforced by Lagrange multipliers. Hence all subdomain stiffness matrices are singular with 
a-priori known kernels which is a great advantage in the numerical solution. With known kernel basis we 
can regularize effectively the stiffness matrix without extra fill in and use any standard sparse Cholesky 
type decomposition method for nonsingular matrices [BDKKM10] . An alternative linear solver based on 
the Schwarz domain decomposition method with overlapping was introduced by [BG94J. 

The structure of the paper is as follows: In Section 2, we introduce a quasistatic scheme of a solid 
mechanics problem. Within this context, we consider the both elastic and elastoplastic models. The elastic 
model is introduced from the methodical purposes because it is an essential part of the investigated elasto- 
plastic model. After its time discretization we summarize the resulting one time step problem. In Section 
3, the finite element space discretization of the one time step problem and its nonlinear algebraic formu- 
lation are described in details. The semismooth Newton method is applied to treat this nonlinearity. The 
TFETI method combined with the projected conjugate gradient algorithm is derived in Section 4 to solve 
the linearized problems appearing in the Newton iterations. Finally, the algorithm for solving the whole 
elastoplastic problem is summarized. In Section 5, the performance of the proposed algorithm is illustrated 
on numerical experiments. Final comments are summarized in Section 6. 

2 Elastic and elastoplastic models 

In this section, we summarize elastic and elastoplastic models in a quasistatic framework. Firstly, we 
introduce basic notation and assumptions that will be used for setting the models. Secondly, we describe 
the linear elastic constitutive model for an isotropic material. Thirdly, we introduce the elastoplastic initial 
value constitutive model based on the described elasticity, the von Mises yield function, the associated 
plastic flow rule and the linear isotropic hardening represented by the accumulated plastic strain. Fourthly, 
we consider time discretization of the constitutive model given by the implicit Euler method. The explicit 
form of the time discretized elastoplastic operator is introduced. 

2.1 Notation and assumptions 

Let us consider a deformable body occupying a domain Q, C M 3 with a Lipschitz continuous boundary 
r = dd. We will describe the state of the body during a loading process by the Cauchy stress tensor 
ff£5, the displacement u G M 3 and the small strain tensor e G S. Here S = M 3 ™ is the space of all 
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symmetric second order tensors. Other variables that are necessary for defining the elastoplastic models will 



be introduced in Subsection 2.3 More details can be found in [NPO08J. 

The above variables depend on the spatial variable and on the time variable t £ Q = [to,t*]. The 

small strain tensor is related to the displacement by the linear relation 



= 2 (V« + (V«) ) • 
The equilibrium equation in the quasistatic case reads 

- div(cr(x, t)) = g(x, t) V(», t) £ £1 x Q, 



(1) 



(2) 



where g(x,t) G M 3 represents the volume force acting at the point x £ Q and the time t G Q. 

Let the boundary V be fixed on a part IV that has a nonzero Lebesque measure with respect to T, i.e. 
we prescribe the homogeneous Dirichlet boundary condition on Tjj: 



u(x,t)=0 V(x,i) G IV x Q. 
On the rest of the boundary IV = T \ IV, we prescribe the Neumann boundary conditions 

a(x, t)n(x) = F(x, t) V(x, t) G T N x Q, 



(3) 



(4) 



where n(x) denotes the exterior unit normal and F(x, t) denotes a prescribed surface forces at the point 
x G IV and the time t G Q. Geometry of £1 with imposed boundary conditions is depicted in Figure [TJ 
Similarly, we can consider other boundary conditions, for example symmetry and periodic conditions. 
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Figure 1: Geometry of the domain Q with imposed boundary conditions. 



For a weak formulation of the investigated problems, it is sufficient to introduce the space of kinematically 
admissible displacements, 

V = {v G [H 1 ^)} 3 : v = on Tu} . (5) 
Then the conditions (^-Q can be written in a weak sense by 

[ (a,e(v)) F dx = [ g T vdx+ [ F T vds G V, Vt G Q. (6) 

Here e(v) is defined by Q, (., .)p and \\-\\f denote the Frobenius scalar product and the corresponding norm 
on the space S, respectively. We assume that the functions a, F, g are sufficiently smooth to be the integrals 
in ^ correctly defined in the Lebesque sense. 

To complete the investigated (generally quasistatic) models, we will prescribe the constitutive relations 
among the stress, the strain and eventually other variables. 
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2.2 Elastic model 

We consider the elastic constitutive model given by the Hooke law for isotropic material, 

a = Ce = \tr(e)I + 2^e (7) 

with the Lame coefficients A, fi. For the sake of simplicity, we assume a homogeneous material, i.e. the 
constant coefficients X, fi > 0. The trace operator of a tensor is denoted by tr(.) and I denotes the identity. 
It will be useful to introduce the volumetric and deviatoric parts of a tensor rj G S by 

vol(ry) = -tr(r/)I, dev(r/) = i] — vol(r/). (8) 
o 

It holds that 

(vol(r ? ),dev(C))F = 0, (devfa),£) F = (devfa),dev(£)>F Vt/,£<e5. (9) 
By ([9]), we can find that the fourth order tensor C, defined by ([7]), is symmetric and elliptic, i.e. 

(C V ,0f = (vX0f, (Cr,,r)) F >2»\\rif F Vr?,£ G 5. (10) 

If we substitute ([7]) into ([6]) , we obtain for any fixed t G Q the weak formulation of the elastic problem. 
Problem 1 (Elastic problem). Find u = u(x,t) G V such that 

a e (u,v)= / g T vdx+ \ F T vds \/v G V, (11) 

where the bilinear form on V reads 

a e (w,v)= / (Ce(w),£(v))pdx 1 w,v£V. (12) 
Jn 

Note that a e (w, v) is symmetric and ^-elliptic by ( |10[ ) and the Korn inequality [NH81j . i.e. 

a e (u;, u) = a e (v, w), 3c > : a e (u, u) > c||u||y Vv,w&V. (13) 



The mentioned properties of the form a e insure that the elastic problem (11) has a unique solution u £V, 
for example by Lax-Milgram lemma |NH81j . Notice that the problem (11) does not depend on the load 
history, so it is a static problem. 

2.3 Elastoplastic initial value constitutive model 

In comparison to elasticity, elastoplasticity is time dependent model where the history of loading is taken 
into account. We will assume associated elastoplasticity with von Mises plastic criterion and linear isotropic 
hardening law (see e.g. |HR99|, IB1991 INPD 08 ) . The elastoplastic initial- value constitutive model consists of 
the following components: 

1. Additive decomposition of the strain tensor into the elastic and plastic parts: 

e = e e + e p . (14) 

2. Linear elastic law between the stress and the elastic strain: 

a = Ce e , (15) 

where the fourth order tensor C is defined by 0. 

3. The von Mises yield function coupled with an isotropic hardening variable k: 

$(g,k) = t/|||dev(<7)|| F - (o-y + H m K) < 0, (16) 
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where a y ,H m > denote the initial yield stress and the hardening modulus, respectively. 

4. The associated plastic flow rule: 

._ . 0$ . [3 dev(a) . ^ ft 

where e v and 7 denote the time derivative of the plastic strain and the plastic multiplier, respectively. 

5. The hardening law based on the accumulated plastic strain rate: 

\¥ p \\f = i. (is) 



Notice that the second equality in (18) follows from (17) 



6. The loading/unloading conditions: 

7>0, $(ct,k)<0, 7 $O,k)=0. (19) 

7. The initial conditions: 

e(x, t ) = £ e (x, t ) = e p (x, t ) = a(x, t ) = 0, k(x, t ) = 0, x £ CI. (20) 

The weak formulation of the corresponding elastoplastic problem can be found in |HR99j . Here we will 
only consider a time discretized elastoplastic model. 

2.4 Time discretized elastoplastic model 

Let us consider the following discretization of the time interval 

t <h < ... <t k < ... <tN = t*. 

Let us denote o k = cr k (x) = a(x,t k ), x G O and similarly for other variables. To approximate the time 
derivatives, we use the implicit Euler method. However, it is also possible to use for example the Crank- 



Nicholson scheme. Then by (14) and (15) 



£ fc+l ~~ £ \ ^ 1 ( (T fc+l ~~ a k+l) 
£ P {tk+l) ~ — -7T, = JT. ) ^k+1 = tk+1 - tk, (21) 

where 

°fc+i = a k + CAe k+ i, Ae k+1 = e k+1 - e k . (22) 



By (14)-(22), we can formulate the time discretized elastoplastic constitutive problem as follows. Given 
the values a k , K k , e k of the stress, the isotropic hardening and the strain, respectively, at the time t k and 
given the incremental strain Ae^+i for the interval [tfc, tfe+i] , solve the following system of algebraic equations 



C (a k+1 -a k+1 ) = A 7fe+ ^- ||devK+i)||F (23) 

K k+ i - Kk = A7 fc+ i (24) 

for the unknowns cr k +i, Kfe+i, and A*y k+ i, subject to the constraints 

A7 fe+ i > 0, ^(fffc+i.Kfc+i) < 0, Aj k+1 <f>{a k+1 ,K k+1 ) = 0. (25) 

This constitutive problem can be solved explicitly by the return mapping concept (see e.g. [B1991 [NPO08J). 
It means that we firstly apply the elastic predictor, i.e. we verify whether $(cr^ +1 , K k ) < 0. If it holds then 

A7 fc+ i = 0, a k+1 = a k+1 . (26) 

If ^{a l k+1 , K k ) > 0, then by the plastic corrector we have 
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where 



re r 



dev(r) 
l|dev(r)|| F ' 



T £ S. 



(28) 



Notice that the second formula in (27) is correctly defined since the denominator ||dev(a^, +1 )||i? > for 
K k) > 0. Let us define the stress and hardening operators T a (T,uj; .) : S — >■ S, T k (t,co; .) : S S 
with respect to parameters t £ S, d £ M+, such that for ry £ S 

3/i [2 



T a (r,u;ri) :-- 
T K (T,uj;rj) :-- 



Cry 



3(i + H„ 



<S> + (t + Ct],oj)h(T + Cry), 



1 



3(i + H, 



-$ + (r + Cr/,w) 



(29) 
(30) 



respectively, where <E> + denotes the positive part of the function <£. Then by (22), (24), (26), (27), (29) and 

(HJ), 

AK fc+ i = T K (a k , «k; Aek+i), Acr fc+ i = T a (a k , K k ; Ae k+ i). (31) 

For the sake of brevity, we will denote the stress operator T a (cr k , n k ; ■) with respect to the current parameters 
a k and K k by T k {.). By [B199, GVOQl [Sy09| , the operator T k : S 1 — > S is potential, Lipschitz continuous, 
strongly monotone, and strongly semismooth on S. 

Let us note that semismoothness was originally introduced by Mifflin |Mi77| for functionals. Qi and J. 
Sun |QS93| extended the definition of semismoothness to vector- valued function to investigate the superlinear 
convergence of the Newton method. The strong semismoothness of the Lipschitz continuous function T k (.) 
means that T k {.) is directionally differentiable on S and has a quadratic approximate property at any rj £ S, 
i.e. for any £ £ S, £ — > 0, and any T£ (ry + £) £ dT k (rj + £), 



T fc (ry + - TMr?) - T fe °(ry + = 0(U\\f)- 



(32) 



Here dT k (r( + £) denotes the set of the Clark generalized derivatives of T k at ry + £. Here we will choose the 
Clark generalized derivative T£ of T k in the following way: 
1. If $(a k + Cry, Kk) < 0, then 

T° k (v) = C 



(33) 



2. If $(a fc + 



, Kk) > 0, then 
T fc °(r?) = C 



3/i 



3/i + H, 
3fi 



i2d$(<j k + Cri,Kk, 
3 drj 



n(a k + Cry) 



where 



3(i + H v 

d$(a k + Ct?, 
drj 

dh(a k + Cry) 



9ry 



2., L „ <9n(<7 fc +Cry) 
-$(<r fc + Cry, — . 



2/uW -n((Tjfc + Cry), 

Id - re(<r fc + Cry) g re(o- fc + Cry) 
M ||dev(a fc +Cr/)|| F 
dev(0, V£ e S. 



(34) 



Notice that T& is not differentiable at ry £ S, &(a k + Cry, K k ) = 0. Otherwise T^(ry) = dTf.[rj)/dr]. 

Let us recall that the stress, strain, hardening and displacement variables also depend on a spatial 
variable x £Vt. We consider the dependence of T k {Ae k ) on x in the following sense: 



T k (Ae k ) = T k {As k ){x) := T a (a k (x), K k (x); Ae k (x)). 



(35) 



Then we can substitute the stress operator T k , defined by (29), into the balance equation ^ to obtain the 
time discretized elastoplastic problem in the incremental form. 
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Problem 2 (One time step elastoplastic problem in the incremental form). Given the stress field a k £ 
[L 2 (Q)]sym and the isotropic hardening field K k £ L 2 (Q,) at the time tk, find the displacement u k+ \ = 
u k + Aujfc+i £ V , where the increment Au k+ \ £ V solves the variational equation 

[ {T k (e(Au k+1 )),e(v)) F dx= [ Agl +1 vdx+ [ AF^vds \/v £ V, (36) 
Jn Jn Jr N 

with loading increments AF k+ i = F k +i — F k , Ag k+ \ = gk+i — 9k- Set the stress and isotropic hardening 
fields (Jfc+i = a k + Acjfc+i, Kk+i = Kk + An k+ i in the next time step tk+\ from the relations 

A<Tfc + i = T a (ak,K k ;e(Au k+ i)), AK k+ i = T K (a k , K k ; e(Au k +i)), (37) 

almost everywhere in Q. 

Problem [2] can be equivalently formulated as a minimization problem |GV091 |Sy09| . Since the operator 



T k is strongly monotone and Lipschitz continuous on S, the non-linear equation ( |36| ) has a unique solution 
Av,k+i £ V (see e.g. [FK80J). As we will see in the next section, we will solve a linearized problem in each 
Newton iteration. To do this, it will be useful to define the bilinear form ak{u) : V x V — > M for u £ V by 

a k (u)(w,v) = / {Tg{e(u))e(w),e(v))dx, v,w£V, (38) 



where the operator Tj!(.) = T£(.)(x) is defined by (33) and (34). Since the operator T k is potential, Lipschitz 
continuous and strongly monotone on S, the bilinear form a k is symmetric and ^-elliptic on V. 



3 Semismooth Newton method in elastoplasticity 

In this section, firstly, we approximate the time discretized elastoplastic problem by the finite element method 
and introduce the corresponding algebraic notation. Secondly, we introduce the semismooth Newton method 
for the problem. 

3.1 Finite element discretization and algebraic formulation 

Details to finite element implementation of elastoplastic problems can be found in |CK02j and |B199| IGV09] . 

For the sake of simplicity, we assume a polynomial 3D domain £7 and use the linear simplex elements. 
The corresponding shape regular triangulation is denoted by Th- Thus the space V is approximated by its 
subspace Vh of piecewise linear and continuous functions. Therefore the spaces of the strains, the stress and 
the isotropic hardening are approximated by piecewise constant functions. 



Similarly to (36), (37), we can formulate the one time step elastoplastic problem after the space dis- 
cretization. Let a kt h, i^k,hi u k,h be piecewise constant stress and hardening variables with respect to the 
triangulation Th at the time t k obtained from a previous time process. 

Problem 3 (One time step elastoplastic problem in the incremental form after space discretization). Find 
the displacement u k +ih = u k,h + ^ u k+i,h £ Vh, where the increment Au k +\h £ Vh solves the variational 
equation 

I (T kjh (£(Au k+ i th )),E(v h )) F dx = / Agl +1 v h dx+ / AF k+l v h ds Vv h eV h , (39) 
Jn Jn Jr N 

where T k ^{.) '■= T a (a kj h, Kk,h> •) • Set the stress and isotropic hardening fields (Jk+i,h = a k,h + ^^k+l,h> 
K k+i,h = Kk t h + Ank+i.h in the next time step t k+ \ from the relations 

Aer fe+lj/l = T a (a k}h , K k ,h; e(A«fc + i^)), A/c fc+lj /j = T K (a kih , Kk,h; e(A«fe + i,/i)) ( 40 ) 

for every elements of Th ■ 

For the sake of simplicity, we do not consider finite element approximation of F k and g k . Similarly 
as in (33) and (34), we can define the generalized derivative T£ h of T k) h and consequently also define the 



approximated bilmear form a k ^h{uh) for Uh £ Vh by 

ak,h( u h)( w h,Vh) = (T k h( £ ( u h))£(w h ),e(v h ))dx, v h ,w h £V h - (41) 



o 
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Each function Vh = {vh,i, Vh,2, Vh,3) £ Vh can be represented by a vector 

v G R", u := Kj(j;i))i 6 {i,„.,^,j 6 {i 1 2,3} ) 

where A/" denotes the number of vertices of the triangulation Th and n = 2>N . The homogeneous Dirichlet 
boundary condition is represented by a restriction matrix Bjj G M mxn , i.e. 

Bt/it = o. (42) 

Let -Rt G M 12xn be a restriction operator for a displacement vector u G M n on a local element T G 7/j, i.e. 

ut = Rtu. (43) 

We denote by u k and Ait fe+1 the displacement vector and the searching displacement increment at the time 
step k, respectively. We denote the load vector represented the volume and surface forces F k ,g k by f k and 
its increment A/ fc . 

Further, we use a vector representation in M 6 of the stress and strain tensors that is typical for an 
implementation of elastic problem, i.e. 



<y = (0"ll,0"22,0'33 ) Cri2,Cr23 ) 0'l3) 



(eiii £22, £33) 2ei2, 2e23> 2ei3) T . 



(44) 



Notice that the stress and strain vectors have different structures in comparison to the above tensor notation. 
Therefore we must carefully distinguish this difference in algebraic representation of the operators T a , T K , 
T kj h, and T kfl . The vectors in sense of stress variables will be denoted by letters er, r, the vectors in sense of 
strain variables will be denoted by letters e, £ p , rj, and £. Let (t^^t an d Kk,T be the algebraic representation 
of Ok,h and Kk,h on an element T £ Th, respectively. 

We introduce the algebraic representations C G M 6x6 , E £ G M 6x6 , E a G M 6x6 , *, n, T K>kiT , T k>T , 
and T° k T of the Hooke tensor C, the deviatoric operator 1^ related to the strain and stress variables, the 
Frobenius norm with respect to a stress variable, the functions <3?, n, and the restrictions of the functions 
T K (cr k \T, Kfclrj •)> T kt h, T kh on T G Th, respectively, with respect to the vector form (44) of the stress and 
strain variables. The forms of matrices C, E £ , E a , and the norm | 
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P:=dm<7(l,l,l,2,2,2), 



and 



r G 



respectively. Consequently the functions 3>, n, T Kjkj T, T k: j>, T kT are defined by (16), (28), (30), (29), (33), 

^3, 



(34) in the following way: 



n(r) 



#(t,k) : = 
E a (r 



\\Eo{T 



1 



3/i + H r 



\E a T\\ a - (oy + H m n) , 



nk,r(v) = n (<r kjT + Crj) 



-* + (<T fc)T + C77, K k T ) 
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3 ix / 2 

T k ,T(v) ■= Cr) - - — —^-\ -* + {(r k ,T + Cr), K k ,T) n k , T (v), 
3/i + H m V 3 

C, if (&k,T + Crj, Kk t j<) < 0, otherwise 
respectively. 

We also introduce the matrix G G M 6x12 representing the algebraical relation between the strain and 
the displacement (the exact form of G is in [ACFK02J), i.e. the strain et on an element T G Th can be 
found by (43) in the form 

e T = GR T u. (45) 
Based on the introduced notation we define the non-linear operator F k : W 1 — > M. n , 

F k (v)=Y J (T k)T {GR T v)) T GR T , v 6 M. n , (46) 
TeT h 



which represents the left hand side in (|39j), and the stiffness matrix K k (v) G M nxn , v G M n , 

K k (v)= £ {Tl T (GR T v)GR T ) T GR T , (47) 
TeT h 

which represents the bilinear form a k) h{vh)- In particular, we denote the matrix K k (ZS.u k+li ) briefly by 
K kj i, where the reason of Z\u k+l i G K™ will be explained in the next section. 
Let 

V := {v G R n \Buv = o} . 



Then by using (46), we can rewrite the equation (39) as follows: find Att fe+1 G V such that 

v T (F k (Au k+1 ) - A/ fc+1 ) = VvtV, (48) 

where A/ fc+1 is the increment of the load vector. Let u k G M. n ~ m , f k G M n_m , K M G rC"-"*) * («-«0 , and 
i^fc : M n ~ m —7- M n ~ m denote the restrictions of u k , f k , K k ^, and F k given by omitting the entries (degrees of 
freedom) corresponding to the prescribed Dirichlet boundary conditions. Then we can rewrite the equation 
( |48[ ) to the following system of non-linear equations: 

find Aii fc+1 G V : F k (Au k+1 ) = Af k+1 . (49) 

The discretized elastoplastic problem can be solved by the following algorithm: 
Algorithm 1 (Solution of discretized elastoplastic problem). 

1: initial step: Uq = o, £o,T = o, <to,t = o, ko,t = o for any T G 7^ 

2: for fc = 0, . . . , N - 1 do 

3: find Aw fc+1 G V: F fc (Aw fc+1 ) = A/ fe+1 

4: for all T G T h do 

5: ^ e fe+l,T = GRxAu k+1 , £ k +l,T = £k,T + A£ i+1T 

6: Act = T fc)T ( Ae fc+1T ), cr fe+1)T = <T fcjT + Act fc+liT 
7: Afc fe+1T = T W)fc)T (Ae fc+ljT ), K fc +i,r = ^fc,T + Ak h1)T 
8: end for 
9: end for 
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3.2 Semismooth Newton method for one time step problem 



The non-linear system of equations (49) is solved by the semismooth Newton method (see e.g. QS93 ). The 
corresponding algorithm is following: 

Algorithm 2 (Semismooth Newton method). 

1: initialization: A-u fc = o 

2: for i = 0,1,2, ... do 

3: find Sui G V: K k)i 8ui = A/ fe+1 - F k (Au ki ) 
4: compute /\u ki+1 = £±u k i + Sui 

5: if l|Au M+1 - Au^H/fllAti^+ill + ||Au fe)i ||) < e Newton then stop 
6: end for 

7: set An H1 = Au fc)i+1 

Here eNewton > is the relative stopping tolerance and Siii G M n_rra is the restriction of Su^ given 
by omitting the entries (degrees of freedom) corresponding to the prescribed Dirichlet boundary conditions. 
The systems of linear equations, which are considered in each Newton iteration, will be solved by the TFETI 
method introduced in the next section. 

In [B199|, IGV091 Sy09| , superlinear local convergence of the algorithm has been derived. Let us note that 



the convergence depends on the discretization parameter h of the triangulation. Therefore we can expect 
that the finer mesh the bigger number of the Newton iterations. In |Sy09|, a damped semismoooth Newton 



method for such a problem has also been described. Such a method has again superlinear local convergence 
and additionally global convergence. 

4 TFETI method for solving linearized problems 

In this section, we describe the TFETI domain decomposition method for solving linearized problems ap- 
pearing in Newton iterations and its optimal solver based on projected conjugate gradient method with a 
lumped preconditioner. Finally, we summarize TFETI based algorithm for solving the whole elastoplastic 
problem. 

4.1 TFETI domain decomposition method 

In the previous section, we showed that we need to solve the following system of the linear equations: 

find Sui G V : K kji 6ui = A/ fe+1 - F k (Au k>i ) (50) 

in each time step k and in each Newton iteration i. Since the TFETI domain decomposition method can be 
described independently of the indices k, i, we will schematically write the problem in the form: 

find u £ V : Ku = /, (51) 

where K, u, f are the restriction of K, u, f with respect to the Dirichlet boundary conditions respectively. 
Let us note that K is a symmetric and positive semidefinite matrix and K is a symmetric and positive 



definite matrix. Therefore the linearized problem (51) has a unique solution. The problem (51) can be 
equivalently rewritten as a minimum problem: 

find u G V : J(u) < J(v), Vi? G V, (52) 

where 

J(v) = ^v T Kv - f T v, veV. 



The corresponding functional representation of (51) is following: find Uh G Vh such that for any Vh G Vh, 



(B h e(u h ),s(v h )) F dx = I g T v h dx+ / F T v h ds- I (T h ,e(v h )) F dx, (53) 



n 
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where restricted on T G Th is a constant fourth order symmetric and elliptic tensor and Th restricted on 
T G 7h is a constant second order tensor from 5". We can see that (53) is related to (50), if we replace 



°h ~ T^ h (Au kjh ,i), r h ~ T k>h (Auk,h,i), F ~ AF k+1 , g ~ Ag k +i, u h ~ 5u h<i 
and set Aitk,h,i and foi^i as the functional representation of Aukj and 5uj, respectively. 



Notice that the problem (53) has a similar scheme as the elastic problem defined in Subsection 2.2 



Therefore the below introduced TFETI method can also be explained for elasticity in the same way. 



The variational problem of the type (53) can be equivalently formulated as a minimization problem: 



find 



Uh G Vh ■ Jh{u h ) < J h {v h ) Mv h G Vh, 



where 



Jh{vh) = \ [ (Ph£(v h ),e(v h )) F dx - 

1 Jn Jn 



g T v h dx 



F v h ds+ / (T h ,e(v h )) F dx. 
r N Jn 



(54) 



(55) 



To apply the TFETI domain decomposition, we tear the body from the part of the boundary with the 
Dirichlet boundary condition, decompose it into subdomains, assign each subdomain by a unique number, 
and introduce new "gluing" conditions on the artificial intersubdomain boundaries and on the boundaries 
with imposed Dirichlet condition (see Figure [2]). 

In particular, the polynomial domain SI is decomposed into a system of s disjoint polynomial subdomains 
f2 p C M 3 , p = 1,2, ... ,s. We assume that the decomposition is consistent with the triangulation Th, i.e. 





a 1 



a 3 



V 



A 



Q 4 



Figure 2: TFETI domain decomposition with subdomain renumbering 



VT G T h 31 p£ {l,2,...,s} : Tcff (56) 

and define 

Tl := {T G T h : T C Cl p }, T h = |J T?. (57) 

P e{i,2,...,s} 

After the decomposition each boundary T p of fP consists of three disjoint parts T^, T P N , and T P G , T p = 
r^ur^ur^-, where 

r p u = r u nv p , r; = r w nF, r p G = |J rg, 

qe{l,2,...,s}\{p} 

with being the part of T p which is glued to £l q , p ^ q. 

Similarly to the definition of Vh, we can define the spaces V p , p = 1, 2, ... s, of piecewise linear and 
continuous approximations of H 1 (fl p ): 

V p := K € HHn p ) : v p h \ T G [Pi] 3 VT G 7?, <| r - = 0}. (58) 
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Let V/j := x Vfi x . . . x V£ be a product space and 

K h :={v h = {vl...,v s h ) £V h : v p h = v q h onT p J Vp, q G {1, 2, . . . , s}, p^q}. 



(59) 



Let us note that we slightly distinguish the notation introduced for the TFETI method from the notation 
introduced in Sections [2] and [3j For example, we write V/j for the TFETI method while in Section [2} we used 
the notation Vh- A similar distinction is also introduced for the below algebraic description of the TFETI 
method. 



Let 



Jh(vh) 



(O h e(vl),e(v p h )) F dx- / g T v p dx 



9.i' 



[ F T v p ds + [ (T h ,e(v p )) F dx 



(60) 



be a functional defined on V^. Then the minimization problem (54) can be equivalently rewritten into the 
form: 

find n h G K h : 3 h (^h) < Jh(v h ) \/v h G K h , (61) 

where u h = ((u h )\ n i,. . . , (u h )\as) and u h G V h solves Q. 

Each function = (y^, v^, . . . , vfy , G V^, can be represented by a vector v G M n , v = u^, . . . , , 
where i> p G IR ra f , ,p G {1,2,..., s}, is the algebraic representation of v p h and n = J2p=i n y Similarly we can 

find the vector f el", f = (fj, /£,..., / J) T , / p G R n p,p G {1, 2, . . . , s}, such that f p is the algebraic 
representation of the load restricted on f2 p and T^. Let the matrix Bg G M mG,xn represent the gluing con- 
ditions introduced in (59) and B[/ G ]R mt7Xn the Dirichlet boundary conditions introduced in (58). More 
details about matrices B G and B[/ can be found in [KVMHDHKC12J. Both matrices can be combined into 
one constraint matrix 



B 



B, 



B G 



m 



Typically m is much smaller than n. Let the matrix K G I 
symmetric positive semidefinite block diagonal matrix, where 



K p = \ T \ {D T GR P T 



nXn , K = diag(K 1 ,K 2 , 



, K s ) denotes a 



Here Dt G IR 6x6 is the algebraic representation of I\|t and Rj, G R 12xn p is a restriction operator for a 
displacement vector u p G M" p to a local element T G 7/f . The diagonal blocks K p , p G {1, 2, . . . , s}, which 
correspond to the subdomains £l p , are positive semidefinite sparse matrices with known kernels, the rigid 
body modes. 



The algebraical formulation of (61) is following: 

find u G V 
J(v) 
V 



J(u) < J(v) Vv G V, 
1 f v T Kv - f T v, 

Bv = o} . 



2 

{vG 



(62) 



Even though (62) is a standard convex quadratic programming problem, its formulation is not suitable 
for numerical solution. The reasons are that K is typically ill-conditioned, singular, and very large. 

The complications mentioned above may be essentially reduced by applying the duality theory of convex 
programming (see, e.g., Dostal |D09j), where all the constraints are enforced by the Lagrange multipliers A. 
The Lagrangian associated with problem ( 62 ) is 



L(v,A) = J(v) + A T Bv. 



It is well known [D09J that (62) is equivalent to the saddle point problem: 

find (u, A) £ 1° x 1' : L(u, v) < L(u, A) < L(v, A) V(v, v) G 



(63) 



(64) 



in sense that u solves (62) if and only if (u, A) solves (64). 
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4.2 Optimal solvers to equality constrained problems 



The solution of (64) leads to the equivalent problem to find (u, A) £ R n x 



satisfying: 



(65) 



The system (65) is uniquely solvable which is guaranteed by the following necessary and sufficient condi- 
tions |BUL05| : 



KerB T = o, 
KerK n KerB 



(66) 
(67) 



n — rank(K) 



Notice that (66) is the condition on the full row-rank of B. Let us mention that an orthonormal basis of 
KerK. is known a-priori and that its vectors are columns of R G M nx ', I 

f — B T A G ImK 



The first equation in ( 65 ) is satisfied if 



(68) 



and 



u = Kt(f-B T A) + Ra (69) 

for an appropriate ael' and arbitrary matrix satisfying KK^K = K. Here is a generalized inverse 
matrix whose application on a vector can be efficiently implemented (see Remark [T]). 

Remark 1. The diagonal block K p , p G {1, 2, . . . , s}, which corresponds to the subdomain is positive 
semidefinite sparse matrix with known kernel basis created by the rigid body modes. This is a great advan- 
tage because all blocks can be effectively regularized without extra fill in and then decomposed using any 
standard sparse Cholesky type factorization method for nonsingular matrices [BDKKM10J. We completely 
avoid problems with zero pivots. The action of on a vector is naturally parallelized with respect to the 
subdomains and computed using backward and forward substitutions. 

Remark 2. Notice that the pseudoinverse can also be applied on vectors which do not belong to 7m(K). 
Therefore we can write K f (f - B T A) = K f f - K t (B T A) in the other text. The random error corresponding 
with this partition is neglicted with respect to the implementation of introduced in Remark [TJ 



The condition (68) can be equivalently written as 



R T (f-B T A) = o. 



Further substituting ( 69 ) into the second equation in ( 65 ) we arrive at 

- BK t B T A + BRa = -BK f f . 

x R l satisfies: 



Summarizing (71) and (70) we find that the pair (A, a) G 



F N T 
N 



A 

a 



(70) 



(71) 



(72) 



where F := BK^B T , N := -R T B T , d := BK+f, and e := -R T f 

Since N is of full row-rank as follow from (67), the inverse (NN 7 )^ 1 exists and Pn := I — N T (NN T ) _1 N 
is well defined and represents the orthogonal projector onto iTerN. Applying Pn on the first equation in 
(72) and checking that PnN t q = o we eliminate a and obtain that A satisfies: 



P N FA = P N d, NA = e. (73) 
A/ m + Xkct into two orthogonal components Xj m G 



In practical computations, we further decompose A 
Im~N T and \xer £ Ker~N, substitute them in to (73) and get the problem 



P N FA Xer = P N (d - FA /m ) on KerN, (74) 

with X Im = N T (NN T )~ e. Equat ion ( |74[ ) is solved efficiently by the projected conjugate gradient method 
with preconditioning (PCGP) [FMR94J using the lumped preconditioner to F in the form F 1 = BKB . 



For more information see |FMR94j . We obtain the following algorithmic scheme for the solution of ( 62 ) : 
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Algorithm 3 (Linear solver based on the TFETI method). 
1: Set N := -R T B T , H := (NN^ -1 , d := BK+f, and e := -R T f. 



2: Compute X] m := N He. 
3: Set d := d - FA/ m . 



4: Compute X Ke r from (74) by PCGP: 



5: r° = d, X° Ker = o. 



6: for j = 1, 2, ... , until convergence do 
7: Project w^ 1 = P^r^ - ■ 




17: Xxer — ^Ker' 



18: Set A := Xj m + A^ e r- 

19: Compute a := HN(d - FA). 

20: Compute u := Kt(f - B T A) + Ret. 

Remark 3. Action of H on a vector may be efficiently implemented by the sparse Cholesky factorization of 



4.3 TFETI based algorithms for solving elastoplastic problem 

In this subsection, we summarize the above proposed algorithms for solving elastoplastic problems and 
modify them with respect to the use of the TFETI domain decomposition method. 

Algorithm 4 (Algorithm for solving elastoplastic problem - sequential version). 

1: uo = o, <t 0i t = o, ko,t = o, T £ Th (initial step) 

2: for k = 0, 1, 2, . . . , N — 1 (time steps) do 

3: set = ° (zero approximation) 

4: for i = 1,2,... (Newton iterations) do 

5: for p = 1,2, ... ,s (cycle over subdomains) do 

6: restrict Aufc+i^-i , Afjt+i into subdomain variables Att^j^j, A/^ +1 

7: call Algorithm [5] with (Auj +li l ,A/j +1 , ffk,T, K>k,T, T £ 7/f) to find output variables K P , i , 

f P k,v Acr fc+i,T, Att fc+ i ]T , T G Tft- 
8: collect K p k i f p k i , into global variables K^, f k>i 

9: end for (cycle over subdomains) 

10: solve by Algorithm [3] the problem: find (5uj 6 V such that 



NN T . 



Jjb,i(*Ui) < Jfc,i(v) Vv e V, with Jfe,j(v) 



1 



v T K M v-f£, 



2 



11: 



Aufc_|_i 5 j = Aufc + i 5 j_i + 5uj (displacement update) 
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12: if || Au fc+M - Aujfc+i^iH/ (||Aujfe +1)i || + HAufc+i^ill) < ^Newton then stop 

13: end for(Newton iter.) 

14: U fe+ i = life + Aufc+1,1, 

15: C"fc + i jT = <T fc;T + A<Tfc + x,Tj *«Jfc+l,T = K fc,T + A/ifc+i^, T G 7ft 

16: end for (cycle over time steps) 



The assembling procedure for subdomain data looks as follows. 
Algorithm 5 (Assemble all data corresponding to a subdomain J7 P ). 

1: Input: Au^ ^j, er fcjT , k A)T , T G 7^. 

2: /fc,t = ^/fc+i 

3: l<l = O 

4: for T G 7^ do 

5: compute |T| (volume of the element T) 
6: A<x fc+liT = T k>T (GR p T Au p k+l i _^ 

7: AK fc+ljT = T kAT (GR p T Au p k+l i _^ 

* fh = fli-\T\{^T) T GR p T 

* = + m {n,T (GR^Aul^) GR p T ) T GR p T 
10: end for (cycle over elements) 

11: Output: K p -, f k -, Aa k+hTl AK k+hT , T G 7^ 



Loop over subdomains and all subdomain operations may be implemented in parallel. Parallelization of 
FETI/TFETI is based on distributing matrix portions among processing units. This allows algorithms to 
be almost the same in sequential and parallel versions; only data structure implementation differs. Most of 
computations (subdomain operations) appearing are purely local and therefore parallelizable without any 
data transfers. Each of cores works with the local part associated with its subdomains. Natural effort 
using the massively parallel computers is to maximize the number of subdomains so that the sizes of the 
subdomain stiffness matrices are reduced. This accelerates their factorization and subsequent pseudoinverse 
application which belongs to the most time consuming action. On the other hand, negative effect of that is 
an increase of the null space dimension and the number of Lagrange multipliers. This leads to larger coarse 
problems, i.e., applications of the projector Pn which are scalable only up to a few thousands of cores and 
then the coarse problem solution starts to dominate. For the numerical solution of such large problems we 
recommend to use a hybrid FETI method |KR10| . 



5 Numerical experiments 

The performance of the proposed algorithms is demonstrated on an elastoplastic homogeneous thin plate 
of sizes 20 x 20 x 2 with the circular hole of radius 1 in the center. Its geometry with imposed boundary 
conditions and indicated symmetry planes is depicted in Figure [3j Thus we consider only eighth of the thin 
plate in our computations (see Figure |4j) and impose symmetry conditions on three symmetry surfaces. A 
similar benchmark was solved in [St03l IGV091 ICKMllj . 

The elastoplastic body Q is made of homegeneous isotropic material with the parameters 

E = 206 900, v = 0.29, a y = 450, and H m = 10000, 

where E and v are the Young modulus and the Poisson ratio, respectively, which are related with the Lame 
coefficients A and /U by the following formulas: 

Ev E 

A = -, -, r = 110 743.8, a = = 801 938. 

(1 + i/)(1-2i/) ' P 2(1 + 1/) 



15 




Figure 3: Geometry of the whole body Figure 4: Geometry of eighth of body 



The indicated traction forces with the history of loading taking into account are prescribed by the function 

g (t) = 400sin(2vrt), t G [to,**], *o = 0, ** = -. 

Since the elastoplastic model is rate-independent and any local unloading is not expected with respect to the 
prescribed load history, the results should be independent of the chosen time discretization. Let us consider 
two variants of the equidistant time discretization characterized by the time step At: 

(a) At = 1/4, N = 1, 

(b) At = 1/32, N = 8. 

For the spatial discretization of 0, let us consider five levels of tetrahedral meshes generated by Ansys with 

520, 1623, 9 947, 166 374, and 309 546 

nodes decomposed into subdomains using Metis. An example of such decomposition is depicted in Figure [5] 
and a zoom near the hole in Figure [6] 




Fi gure 5: Domain decomposition into 10 subdomains Figure 6: Zoom of Figure [5] near the hole 

The proposed algorithms were implemented in MatSol library [KMBKVD09J developed in Matlab and 
parallelized using Matlab Distributed Computing Server and Matlab Parallel Toolbox. For all computations 
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we use maximum 28 cores with 2GB memory per core of the HP Blade system, model BLc7000. The 
stopping tolerances of the Newton and the PCGP algorithms are 



e Newton = 10 and epcGP = 10 



-7 



(75) 



respectively. 

In Table [TJ we report the numerical results for the time discretization (a) and different mesh levels. The 
numbers of subdomains are chosen to keep the number of nodes per subdomain approximately constant 
except the coarsest mesh level. We see that the number of the Newton iterations remains almost constant 
for all meshes. Similar behavior was observed in the numerical results by [GV09]. Note that even the 
number of the PCGP iterations increases only moderately for larger meshes. In Table [2j the number of 
the PCGP iterations, the number of plastic elements, and the relative convergence criterion are reported 
for each Newton's iteration of Algorithm [4j the time discretization (a), and the finest mesh level. In this 
case, we set the stopping tolerances e^ewton = 10~ 9 and epccp = 10~ 14 so accurate to observe quadratic 
convergence of the Newton method after identification of the plastic zone. Such behavior agrees with the 
theoretical results in [B199] . 

The computational history of Algorithm [4] for the time discretization (b), the finest mesh level, and the 

3 The corresponding development of the plastic zone is 
7] - [TUt respectively. The growing zone of plastic elements 



stopping tolerances (75) is documented in Table 



depicted for the times t^, t±, t%, and ts in Figures 
results from the monotonically increased loading. Distributions of the von Mises stress ||dev(cr)||_p and the 



total displacement [ u| at the final time t* are depicted in Figures 11 and 12 respectively. 



Mesh level 


1 


2 


3 


4 


5 


Number of mesh nodes 


520 


1 623 


9 947 


166 374 


309 546 


Number of mesh elements 


1 441 


6 279 


48 287 


931 709 


1 758 907 


Number of subdomains 


1 


1 


8 


135 


270 


Number of CPU cores 


2 


2 


9 


28 


28 


Number of primal variables 


1 560 


4 869 


33 414 


623 904 


1 183 101 


Number of dual variables 


284 


618 


5 448 


136 904 


272 053 


Number of plastic elements 


252 


2 152 


20 790 


420 760 


796 520 


Number of Newton iterations. 


5 


6 


6 


6 


6 


Total number of PCGP iterations 


237 


391 


573 


718 


724 


Total time in seconds 


8 


60 


77 


928 


1796 



Table 1: Numerical results of Algorithm [4] for different mesh levels and the time discretization (a) 



Newton iter. 


Number of 


Number of 


Stopping 




PCGP iters. 


plastic elem. 


criterion 


1 


224 





1 


2 


312 


808 116 


1.4016e-l 


3 


310 


797 081 


4.6051e-2 


4 


319 


788 802 


6.3207e-3 


5 


318 


795 245 


4.5604e-4 


6 


319 


796 519 


1.0735e-5 


7 


321 


796 596 


1.1790e-8 


8 


314 


796 596 


9.0488e-14 



Table 2: Computational history of Algorithm S for the time discretization (a) and the finest mesh level 



Remark 4. The numbers of the PCGP iterations may be decreased by reducing the required accuracy or by 
assuming the inexact Newton method [B199J. 

Comparing time discretizations (a) and (6) we see that the resulting number of plastic elements differs 
at the final time t* . The difference is less than 0.5% and is caused by the roundoff errors, the use of iterative 
solvers, and the numerical evaluation of the yield function which decides whether an element plasticizes or 
not. 
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Figure 7: Elastic (gray color) and plastic (black color) 
elements at time t 2 



Figure 8: Elastic (gray color) and plastic (black color) 
elements at time t± 




Figure 9: Elastic (gray color) and plastic (black color) 
elements at time i@ 



- 573.0 




Figure 11: Distribution of von Mises stress ||dev(cr)||f 
at f 8 




Figure 10: Elastic (gray color) and plastic (black color) 
elements at time t§ 



-0.02218 




Figure 12: Total displacement ||u|| at t% 
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1 imc 


Number of 


Total number of 


Number of 


Solution 


step 


Newton iters. 


rbiar itciations 


plastic elems. 


time 


1 


o 
z 


zzo 


U 


oob 


o 
z 


o 
z 


nno 
ZZO 


o /UU 


CQ /I 

0o4 


Q 
O 


4 


4aY 


lol 001 


1 1 KO 

1 loz 


A 

4 


4 


coo 

ozz 


316 718 


1 OO 7 
1 237 


K 

o 






OUZ UUJ 




6 


5 


863 


671 067 


1 678 


7 


5 


899 


762 501 


1 638 


8 


4 


662 


799 989 


1 445 



Table 3: Computational history of Algorithm H for the time discretization (6) and the finest mesh level 



6 Conclusion 

We have proposed the algorithm for the efficient parallel implementation of elastoplastic problems with the 
von Mises plastic criterion and the linear isotropic hardening law which is based on the TFETI domain 
decomposition method. For the time discretization we used implicit Euler method and for the space dis- 
cretization of the respective one time step elastoplastic problem the finite element method. The latter results 
in a system of nonlinear equations with a strongly semismooth and strongly monotone operator. Thus the 
semismooth Newton method was applied and respective linearized problems were solved in parallel using 
TFETI. 

The performance of our algorithm was demonstrated on the 3D elastoplastic thin plate with the hole 
in the center and prescribed loading history. Numerical results for different time discretizations and mesh 
levels were presented and discussed. Local quadratic convergence of the semismooth Newton method was 
observed after identification of the plastic zone which is in accordance with the theoretical results. Our 
future plan is to apply the proposed TFETI based algorithm on elastoplastic multi-body contact problems 
of mechanics. This has already been done successfully for the pure elastic case |DKVBM101 lDKMBVH12j. 
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